## ---------------------------------------------------------------------- #
## Figure 9
## Christopher J. Fariss
## "Yes, Human Rights Practices Are Improving Over Time"
## American Political Science Review
## ---------------------------------------------------------------------- #
rm(list = ls())

library(plotrix)

VDEM <- read.csv("V-Dem-DS-CY-v7.csv")
VDEM <- subset(VDEM, select=c(COWcode, year, v2cltort, v2cltort_sd, v2clkill, v2clkill_sd, v2x_clphy, v2x_clphy_codehigh, v2x_clphy_codelow, v2csreprss, v2csreprss_sd, country_name))
names(VDEM) <- c("COW", "YEAR", "v2cltort", "v2cltort_sd", "v2clkill", "v2clkill_sd", "v2csreprss", "v2csreprss_sd", "country_name")



year <- 1946:2016
data <- subset(VDEM, YEAR>=1946 & !is.na(v2cltort))

TORTURE <- FALSE
KILL <- FALSE
REPRESS <- TRUE

if(TORTURE) data$latentmean <- data$v2cltort
if(TORTURE)data$latentsd <- data$v2cltort_sd

if(KILL) data$latentmean <- data$v2clkill
if(KILL)data$latentsd <- data$v2clkill_sd

if(REPRESS) data$latentmean <- data$v2csreprss
if(REPRESS) data$latentsd <- data$v2csreprss_sd

temp1 <- temp2 <- NA
for(i in 1:length(year)){
    m <- data$latentmean[data$YEAR==year[i]]
    v <- data$latentsd[data$YEAR==year[i]]
    temp <- apply(matrix(rnorm(matrix(1,ncol=1000, nrow=length(m)), mean=m, sd=v),ncol=1000, nrow=length(m)), 2, mean)
    temp1[i] <- mean(temp)
    temp2[i] <- sd(temp)
}
latentmean <- temp1
latentsd <- temp2
YEAR <- year

data <- as.data.frame(cbind(YEAR, latentmean, latentsd))

COLOR <- c(2,3)
COLOR <- c(1,1)

temp <- matrix(-999, nrow=5000, ncol=length(year))



temp <- matrix(-999, nrow=5000, ncol=length(year))


for(i in year-1945){
    temp[,i]  <- rnorm(5000, mean=data$latentmean[data$YEAR==year[i]], sd=data$latentsd[data$YEAR==year[i]])
}


# calculate the mean value of each column in the matrix; i.e., the mean of each column of 5000 draws
point.estimates <- apply(temp, 2, mean)


par(mar=c(4,5,4,2), cex=1.2, font=1, font.lab=1)

#COLOR <- "navy"
COLOR <- "black"

# plot the country-year point.estimates
plot(point.estimates, ylim=c(-0.6,1.60), yaxt="n", xaxt="n", pch=19, col=grey(.25), xlab="", ylab="", cex=.8, type="n", xlim=c(1, 70))
axis(1, at=(1:70), label=rep("", length(1:70)))
axis(1, at=(seq(5,65,10)), label=seq(1950,2010,10), lwd=1.5)
axis(2, at=c(-0.75,-0.50,-.25,0,0.25,0.50,0.75, 1.0, 1.25, 1.5), las=2)
#axis(2, at=c(-1,-0.5,0,0.5,1.0), las=2)
if(KILL) mtext(side=2, "Mean V-DEM Killing", line=3.5, cex=1.7)
if(TORTURE) mtext(side=2, "Mean V-DEM Torture", line=3.5, cex=1.7)
if(REPRESS) mtext(side=2, "Mean V-DEM CSO Repression", line=3.5, cex=1.7)
mtext(side=4, "More Abuse", font=1, line=0.3, cex=1.7, at=-0.10)
mtext(side=4, "More Respect", font=1, line=0.3, cex=1.7, at=0.90)
arrows(75.25, -.50, 75.25, -.65, xpd = TRUE, lwd=3, length=.15)
arrows(75.25, 1.4, 75.25, 1.6, xpd = TRUE, lwd=3, length=.15)

abline(h=-.50,lty=2, col=grey(.75), lwd=2)
abline(h=-.25,lty=2, col=grey(.75), lwd=2)
abline(h=0,lty=2, col=grey(.75), lwd=2)
abline(h=.25,lty=2, col=grey(.75), lwd=2)
abline(h=.5,lty=2, col=grey(.75), lwd=2)
abline(h=.75,lty=2, col=grey(.75), lwd=2)
abline(h=1.00,lty=2, col=grey(.75), lwd=2)
abline(h=1.25,lty=2, col=grey(.75), lwd=2)


for(i in 1:length(year)){
    lines(c(i,i), c(apply(temp, 2, quantile, .025)[i], apply(temp, 2, quantile, .975)[i]), lwd=2, col=grey(.6))
}

points(point.estimates, col=COLOR, pch=21, bg=grey(.75), cex=.8)

if(KILL) mtext(side=3, "Estimated Yearly Average of V-DEM \nState Sanctioned Killing Variable", line=0.3, cex=1.7)
if(TORTURE) mtext(side=3, "Estimated Yearly Average of V-DEM \nState Sanctioned Torture Variable", line=0.3, cex=1.7)
if(REPRESS) mtext(side=3, "Estimated Yearly Average of V-DEM \nState Sanctioned CSO Repression Variable", line=0.3, cex=1.7)
